# Following code provided by Preacher and Selig (2012)

a1 =  0.226 # estimated coefficient JS on openness
a2 =  0.192 # estimated coefficient JS on conscientiousness
a3 =  0.0384 # estimated coefficient JS on extraversion
a4 =  0.104 # estimated coefficient JS on agreeableness
a5 =  0.00335 # estimated coefficient JS on neuroticism

a1std = 0.0935 # SE of coefficient a1
a2std = 0.0833 # SE of coefficient a2
a3std = 0.0812 # SE of coefficient a3
a4std = 0.0702 # SE of coefficient a4
a5std = 0.0627 # SE of coefficient a5
  
third.b1 =  0.0214 # estimated coefficient openness on best in third grade
third.b2 =  0.0716 # estimated coefficient conscientiousness on best in third grade
third.b3 =  -0.0431 # estimated coefficient extraversion on best in third grade
third.b4 =  -0.0160 # estimated coefficient agreeableness on best in third grade
third.b5 =  -0.0361 # estimated coefficient neuroticism on best in third grade
  
third.b1std = 0.0278 # SE on coefficient b1
third.b2std = 0.0287 # SE on coefficient b2
third.b3std = 0.0344 # SE on coefficient b3
third.b4std = 0.0286 # SE on coefficient b4
third.b5std = 0.0336 # SE on coefficietn b5

elem.b1 = 0.0291 # estimated coefficient openness on best in elementary
elem.b2 = 0.0282 # estimated coefficient conscientiousness on best in elementary
elem.b3 = 0.0218 # estimated coefficient extraversion on best in elementary
elem.b4 = -0.0266 # estimated coefficient agreeableness on best in elementary
elem.b5 = 0.0658 # estimated coefficient neuroticism on best in elementary

elem.b1std =  0.0323 # SE on coefficient b1
elem.b2std =  0.0249 # SE on coefficient b2
elem.b3std =  0.0288 # SE on coefficient b3
elem.b4std =  0.0322 # SE on coefficient b4
elem.b5std =  0.0311 # SE on coefficient b5

top.b1 =  0.0584 # estimated coefficient openness on probability in top section
top.b2 =  0.00269 # estimated coefficient conscientiousness on probability in top section
top.b3 =  0.0277 # estimated coefficient extraversion on probability in top section
top.b4 =  0.00733 # estimated coefficient on agreeableness on probability in top section
top.b5 =  0.0170 # estimated coefficient on neuroticims on probability in top section

top.b1std = 0.0201 # SE on coefficient b1
top.b2std = 0.0246 # SE on coefficient b2
top.b3std = 0.0256 # SE on coefficient b3
top.b4std = 0.0256 # SE on coefficient b4
top.b5std = 0.0252 # SE on coefficient b5

enroll.b1 = 0.0171 # estimated coefficient openness on current enroll 
enroll.b2 = 0.0283 # estimated coefficient conscientiousness on current enroll
enroll.b3 = 0.0150 # estimated coefficient extraversion on current enroll
enroll.b4 = -0.0177 # estimated coefficient agreeableness on current enroll
enroll.b5 = -0.0327 # estimated coefficient neuriticism on current enroll

enroll.b1std =  0.0184 # SE on coefficient b1
enroll.b2std =  0.0129 # SE on coefficient b2
enroll.b3std =  0.0134 # SE on coefficient b3
enroll.b4std =  0.0162 # SE on coefficient b4 
enroll.b5std =  0.0165 # SE on coefficient b5
  
rep = 20000 # number of simulated values
conf = 95 # confidence level

a1vec = rnorm(rep)*a1std+a1 # create vector of simulated a1 coefficients
a2vec = rnorm(rep)*a2std+a2 # create vector of simulated a2 coefficients
a3vec = rnorm(rep)*a3std+a3 # create vector of simulated a3 coefficients
a4vec = rnorm(rep)*a4std+a4 # create vector of simulated a4 coefficients
a5vec = rnorm(rep)*a5std+a5 # create vector of simulated a5 coefficients

third.b1vec = rnorm(rep)*third.b1std+third.b1 # create vector of simulated b1 coefficients
third.b2vec = rnorm(rep)*third.b2std+third.b2 # create vector of simulated b2 coefficients
third.b3vec = rnorm(rep)*third.b3std+third.b3 # create vector of simulated b3 coefficients 
third.b4vec = rnorm(rep)*third.b4std+third.b4 # create vector of simulated b4 coefficients
third.b5vec = rnorm(rep)*third.b5std+third.b5 # create vector of simulated b5 coefficients

elem.b1vec = rnorm(rep)*elem.b1std+elem.b1
elem.b2vec = rnorm(rep)*elem.b2std+elem.b2 
elem.b3vec = rnorm(rep)*elem.b3std+elem.b3
elem.b4vec = rnorm(rep)*elem.b4std+elem.b4
elem.b5vec = rnorm(rep)*elem.b5std+elem.b5

top.b1vec = rnorm(rep)*top.b1std+top.b1
top.b2vec = rnorm(rep)*top.b2std+top.b2 
top.b3vec = rnorm(rep)*top.b3std+top.b3
top.b4vec = rnorm(rep)*top.b4std+top.b4
top.b5vec = rnorm(rep)*top.b5std+top.b5

enroll.b1vec = rnorm(rep)*enroll.b1std+enroll.b1
enroll.b2vec = rnorm(rep)*enroll.b2std+enroll.b2 
enroll.b3vec = rnorm(rep)*enroll.b3std+enroll.b3
enroll.b4vec = rnorm(rep)*enroll.b4std+enroll.b4
enroll.b5vec = rnorm(rep)*enroll.b5std+enroll.b5

# Openness on best in third grade
openness_third = (a1vec*third.b1vec)

low = (1-conf/100)/2
upp = ((1-conf/100)/2)+(conf/100)
LL = quantile(openness_third,low) # lower limit of confidence interval
UL = quantile(openness_third,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
openness.third.fig <- hist(openness_third,breaks = 'FD',col = 'skyblue', 
     xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
     main = '')

# Conscientiousness on best in third grade
conscientiousness_third = (a2vec*third.b2vec)

LL = quantile(conscientiousness_third,low)
UL = quantile(conscientiousness_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
conscientiousness.third.fig <- hist(conscientiousness_third,breaks = 'FD',col = 'skyblue',
     xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
     main = '')

# Extraversion on best in third grade
extraversion_third = (a3vec*third.b3vec)

LL = quantile(extraversion_third,low)
UL = quantile(extraversion_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
extraversion.third.fig <- hist(extraversion_third,breaks = 'FD',col = 'skyblue',
     xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
     main = '')

# Agreeableness on best in third grade
agreeableness_third = (a4vec*third.b4vec)

LL = quantile(agreeableness_third,low)
UL = quantile(agreeableness_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
agreeableness.third.fig <- hist(agreeableness_third,breaks = 'FD',col = 'skyblue',
     xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
     main = '')

# Neuroticism on best in third grade
neuroticism_third = (a5vec*third.b5vec)

LL = quantile(neuroticism_third,low)
UL = quantile(neuroticism_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
neuroticism.third.fig <- hist(neuroticism_third,breaks = 'FD',col = 'skyblue',
     xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
     main = '')

# Big 5 on best in third grade
total_third = (a1vec*third.b1vec) + (a2vec*third.b2vec) + (a3vec*third.b3vec) + (a4vec*third.b4vec) + (a5vec*third.b5vec)

LL = quantile(total_third,low)
UL = quantile(total_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
total.third.fig <- hist(total_third,breaks = 'FD',col = 'skyblue',
     xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
     main = '')

# Openness on best in elementary
openness_elem = (a1vec*elem.b1vec)

LL = quantile(openness_elem,low) # lower limit of confidence interval
UL = quantile(openness_elem,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
openness.elem.fig <- hist(openness_elem,breaks = 'FD',col = 'skyblue', 
                           xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                           main = '')

# Consientiousness on best in elementary
conscientiousness_elem = (a2vec*elem.b2vec)

LL = quantile(conscientiousness_elem,low) # lower limit of confidence interval
UL = quantile(conscientiousness_elem,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
conscientiousness.elem.fig <- hist(conscientiousness_elem,breaks = 'FD',col = 'skyblue', 
                          xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                          main = '')

# Extraversion on best in elementary
extraversion_elem = (a3vec*elem.b3vec)

LL = quantile(extraversion_elem,low) # lower limit of confidence interval
UL = quantile(extraversion_elem,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
extraversion.elem.fig <- hist(extraversion_elem,breaks = 'FD',col = 'skyblue', 
                          xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                          main = '')

# Agreeableness on best in elementary
agreeableness_elem = (a4vec*elem.b4vec)

LL = quantile(agreeableness_elem,low) # lower limit of confidence interval
UL = quantile(agreeableness_elem,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
agreeableness.elem.fig <- hist(agreeableness_elem,breaks = 'FD',col = 'skyblue', 
                          xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                          main = '')

# Neuroticism on best in elementary
neuroticism_elem = (a5vec*elem.b5vec)

LL = quantile(neuroticism_elem,low) # lower limit of confidence interval
UL = quantile(neuroticism_elem,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
neuroticism.elem.fig <- hist(neuroticism_elem,breaks = 'FD',col = 'skyblue', 
                          xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                          main = '')

# Big 5 on best in elementary
total_elem = (a1vec*elem.b1vec) + (a2vec*elem.b2vec) + (a3vec*elem.b3vec) + (a4vec*elem.b4vec) + (a5vec*elem.b5vec)

LL = quantile(total_elem,low)
UL = quantile(total_elem,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
total.elem.fig <- hist(total_elem,breaks = 'FD',col = 'skyblue',
                        xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                        main = '')



# Openness on probability in top section
openness_top = (a1vec*top.b1vec)

low = (1-conf/100)/2
upp = ((1-conf/100)/2)+(conf/100)
LL = quantile(openness_top,low) # lower limit of confidence interval
UL = quantile(openness_top,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
openness.top.fig <- hist(openness_top,breaks = 'FD',col = 'skyblue', 
                           xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                           main = '')

# Conscientiousness on probability in top section
conscientiousness_top = (a2vec*top.b2vec)

LL = quantile(conscientiousness_top,low)
UL = quantile(conscientiousness_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
conscientiousness.top.fig <- hist(conscientiousness_top,breaks = 'FD',col = 'skyblue',
                                    xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                                    main = '')

# Extraversion on probability in top section
extraversion_top = (a3vec*top.b3vec)

LL = quantile(extraversion_top,low)
UL = quantile(extraversion_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
extraversion.top.fig <- hist(extraversion_top,breaks = 'FD',col = 'skyblue',
                               xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                               main = '')

# Agreeableness on probability in top section
agreeableness_top = (a4vec*top.b4vec)

LL = quantile(agreeableness_top,low)
UL = quantile(agreeableness_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
agreeableness.top.fig <- hist(agreeableness_top,breaks = 'FD',col = 'skyblue',
                                xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                                main = '')

# Neuroticism on probability in top section
neuroticism_top = (a5vec*top.b5vec)

LL = quantile(neuroticism_top,low)
UL = quantile(neuroticism_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
neuroticism.top.fig <- hist(neuroticism_top,breaks = 'FD',col = 'skyblue',
                              xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                              main = '')

# Big 5 on best in third grade
total_top = (a1vec*top.b1vec) + (a2vec*top.b2vec) + (a3vec*top.b3vec) + (a4vec*top.b4vec) + (a5vec*top.b5vec)

LL = quantile(total_top,low)
UL = quantile(total_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
total.top.fig <- hist(total_top,breaks = 'FD',col = 'skyblue',
                        xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                        main = '')



# Openness on current enrollment
openness_enroll = (a1vec*enroll.b1vec)

low = (1-conf/100)/2
upp = ((1-conf/100)/2)+(conf/100)
LL = quantile(openness_enroll,low) # lower limit of confidence interval
UL = quantile(openness_enroll,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
openness.enroll.fig <- hist(openness_enroll,breaks = 'FD',col = 'skyblue', 
                           xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                           main = '')

# Conscientiousness on current enrollment
conscientiousness_enroll = (a2vec*enroll.b2vec)

LL = quantile(conscientiousness_enroll,low)
UL = quantile(conscientiousness_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
conscientiousness.enroll.fig <- hist(conscientiousness_enroll,breaks = 'FD',col = 'skyblue',
                                    xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                                    main = '')

# Extraversion on current enrollment
extraversion_enroll = (a3vec*enroll.b3vec)

LL = quantile(extraversion_enroll,low)
UL = quantile(extraversion_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
extraversion.enroll.fig <- hist(extraversion_enroll,breaks = 'FD',col = 'skyblue',
                               xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                               main = '')

# Agreeableness on current enrollment
agreeableness_enroll = (a4vec*enroll.b4vec)

LL = quantile(agreeableness_enroll,low)
UL = quantile(agreeableness_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
agreeableness.enroll.fig <- hist(agreeableness_enroll,breaks = 'FD',col = 'skyblue',
                                xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                                main = '')

# Neuroticism on current enrollment
neuroticism_enroll = (a5vec*enroll.b5vec)

LL = quantile(neuroticism_enroll,low)
UL = quantile(neuroticism_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
neuroticism.enroll.fig <- hist(neuroticism_enroll,breaks = 'FD',col = 'skyblue',
                              xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                              main = '')

# Big 5 on currently enrollment
total_enroll = (a1vec*enroll.b1vec) + (a2vec*enroll.b2vec) + (a3vec*enroll.b3vec) + (a4vec*enroll.b4vec) + (a5vec*enroll.b5vec)

LL = quantile(total_enroll,low)
UL = quantile(total_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
total.enroll.fig <- hist(total_enroll,breaks = 'FD',col = 'skyblue',
                        xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                        main = '')

